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ABSTRACT 

We have calculated the cosmic ray (CR) acceleration at young remnants from Type la supernovae expanding 
into a uniform interstellar medium (ISM). Adopting quasi-parallel magnetic fields, gasdynamic equations and the 
diffusion convection equation for the particle distribution function are solved in a comoving spherical grid which 
expands with the shock. Bohm-type diffusion due to self-excited Alfven waves, drift and dissipation of these waves 
in the precursor and thermal leakage injection were included. With magnetic fields amplified by the CR streaming 
instability, the particle energy can reach up to 10 1G Z eV at young supernova remnants (SNRs) of several thousand 
years old. The fraction of the explosion energy transferred to the CR component asymptotes to 40-50 % by that 
time. For a typical SNR in a warm ISM, the accelerated CR energy spectrum should exhibit a concave curvature 
with the power-law slope flattening from 2 to 1.6 at E > 0.1 TeV. 
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I. INTRODUCTION 

Diffusive shock acceleration (DSA) is widely ac- 
cepted as the primary mechanism through which cos- 
mic rays are generated in a wide range of astrophysical 
shocks (Drury, 1983; Malkov & Drury 2001 and ref- 
erences therein; Kang & Jones 2002; Kang 2003). It 
is well known that the CR energy density is compara- 
ble to the gas thermal energy density in the ISM and 
plays important dynamical roles in the evolution of our 
Galaxy. Most of galactic cosmic rays, at least up to 
10 14 eV of the proton energy, are believed to be accel- 
erated by SNRs within our Galaxy via Fermi first order 
process (Blanford & Ostriker 1978; Lagage & Cesarsky 
1983; Blandford & Eichler 1987; Drury et al. 2001). 

Time-dependent, kinetic simulations of the CR ac- 
celeration at SNRs have shown that up to 50 % of explo- 
sion energy can be converted to CRs, when a fraction 
10~ 4 — 10~ 3 of incoming thermal particles are injected 
into the CR population at the subshock (e.g., Berezhko, 
Ksenofontov, & Yelshin 1995; Berezhko, & Volk 1997; 
Volk & Berezhko 2005; Kang & Jones 2006). This 
should be enough to replenish the galactic CRs escaping 
from our Galaxy with Lcr ~ 10 41 erg s _1 . X-ray obser- 
vations of young SNRs such as SN1006 and Cas A indi- 
cate the presence of TeV electrons emitting nonthermal 
synchrotron emission immediately inside the outer SNR 
shock (Koyama et al. 1995; Bamba et al. 2003). They 
provide clear evidences for the efficient acceleration of 
the CR electrons at SNR shocks. Also recent HESS 
observation of SNR RXJ1713. 7-3946 indicates possible 
detections of ir° decay 7 rays from the hadronic CRs 
accelerated by the SNR shock (Aharonian et al. 2004; 
Berezhko & Volk 2006). 

In the kinetic equation approach to numerical study 
of CR acceleration at shocks, the diffusion-convection 



equation for the particle momentum distribution, f(p), 
is solved along with suitably modified gasdynamic equa- 
tions (e.g., Kang & Jones 1991). Accurate solutions to 
the CR diffusion-convection equation require a com- 
putational grid spacing significantly smaller than the 
particle diffusion length, Ax <C x<i(p) = K (p)/u s , 
where n(p) is diffusion coefficient and u s is the shock 
speed. In a realistic diffusion transport model, the dif- 
fusion coefficient has a steep momentum dependence, 
k(p) oc p s , with s ~ 1 — 2. So a wide range of length 
scales is required to be resolved in order to follow 
the CR acceleration from the injection energy (typ- 
ically Pinj/nipC ~ 10~ 2 ) to highly relativistic energy 
(p/m p c 3> 1). This constitutes an extremely challeng- 
ing numerical task, requiring rather extensive compu- 
tational resources. 

To overcome this numerical problem, Berezhko and 
collaborators (e.g., Berezhko et al. 1995) introduced a 
"change of variables technique" in which the radial co- 
ordinate is normalized by the diffusion length, Xd(p), 
at each particle momentum for the upstream region. 
This allowed them to solve the coupled system of gas- 
dynamic equations and the CR transport equation with 
k(p) ex p. Their scheme was designed for simulations of 
supernova remnants, which were represented by piston- 
driven spherical shocks in one-dimensional (ID) spher- 
ical geometry. 

On the other hand, Kang and collaborators have 
taken an alternative approach that is based on a more 
conventional Eulerian formalism. Adaptive Mesh Re- 
finement (AMR) technique and subgrid shock tracking 
technique were combined to build CRASH (Cosmic- 
Ray Amr SHock) code in ID plane-parallel geometry 
(Kang et al. 2001) and in ID spherical symmetric ge- 
ometry (Kang & Jones 2006). In order to implement 
the shock tracking and AMR techniques effectively in 
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a spherical geometry, we solve the fluid and diffusion- 
convection equations in a frame comoving with the 
outer spherical shock. Adopting a comoving frame 
turns out to be a great numerical success, since we 
can achieve numerical convergence at a grid resolution 
much coarser than that required in an Eulcrian grid. In 
the comoving grid, the shock remains at the same loca- 
tion, so the compression rate is applied consistently to 
the CR distribution at the subshock, resulting in much 
more accurate and efficient low energy CR acceleration. 

In the present paper, we apply the spherical CRASH 
code for the CR acceleration at remnant shocks from 
Type la SNe expanding into a uniform interstellar 
medium, assuming a quasi-parallel field geometry. De- 
tails of the numerical method are described in §11. The 
simulation results are presented and discussed in §111, 
followed by a summary in §IV. 

II. NUMERICAL METHOD 
(a) BASIC EQUATION 

Here we consider the CR acceleration at a quasi- 
parallel shock where the magnetic field lines are par- 
allel to the shock normal. So we solve the standard 
gasdynamic equations with CR pressure terms added 
in the Eulerian formulation for one dimensional spher- 
ical symmetric geometry. 
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S + ii(Su) = ~Su+ [W(r, t) - L(r, t)} , 



dt dr 



(4) 

where P g and P c are the gas and the CR pressure, 
respectively, e g = Pg/[p(j g — 1)] +u 2 /2 is the total en- 
ergy of the gas per unit mass. The evolution of a mod- 
ified entropy, S = Pg/p 7 " -1 , is followed everywhere 
except across the subshock, since for strongly shocked 
flows numerical errors in computing the gas pressure 
from the total energy can lead to spurious entropy gen- 
eration with standard methods, especially in the shock 
precursor (Kang, Jones, & Gieseler 2002). Total energy 
conservation is applied only across the subshock. The 
remaining variables, except for L and W, have stan- 
dard meanings. The injection energy loss term, L(r, t), 
accounts for the energy carried by the suprathermal 
particles injected into the CR component at the sub- 
shock. Gas heating due to Alfven wave dissipation in 



the upstream region is represented by the term 

BP 



(5) 



where va = B/^Airp is the Alfven speed. This term 
derives from a simple model in which Alfven waves are 
amplified by streaming CRs and dissipated locally as 
heat in the precursor region (e.g., Jones 1993). 

The CR population is evolved by solving the diffusion- 
convection equation, 



dg , 
dt +iu 
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where g — p 4 f, with f(p, r, t) the pitch angle averaged 
CR distribution, and y = ln(p), while n(r,y) is the 
diffusion coefficient parallel to the field lines (Skilling 
1975). For simplicity we express the particle momen- 
tum, p in units m p c hereafter and consider only the 
proton CR component. The wave speed is set to be 
Uw = V A m the upstream region, while we use u w = 
in the downstream region. This term reflects the fact 
that the scattering by Alfven waves tends to isotropizc 
the CR distribution in the wave frame rather than the 
gas frame. 

(b) Spherical CRASH code 

Details of the CRASH code in ID spherical sym- 
metric geometry can be found in Kang & Jones (2006), 
so we briefly describe the basic features here. We solve 
the equations (l)-(6) in a comoving frame that expands 
with the instantaneous shock speed, since a spherical 
shock can be made to be stationary by adopting comov- 
ing variables which factor out a uniform expansion or 
contraction. Because the shock is at rest and tracked 
accurately as a true discontinuity, we can refine the 
region around the gas subshock at an arbitrarily fine 
level. The AMR technique allows us to "zoom in" in- 
side the precursor structure with a hierarchy of small, 
refined grid levels applied around the shock. The re- 
sult is an enormous savings in both computational time 
and data storage over what would be required to solve 
the problem using more traditional methods on a single 
fine grid. 

(c) The Thermal Leakage Injection Model 

The injection rate with which suprathermal particles 
are injected into CRs at the subshock depends in gen- 
eral upon the shock Mach number, field obliquity angle, 
and strength of Alfven turbulence responsible for scat- 
tering. The CRASH codes treat this process naturally 
and self-consistently via "thermal leakage" through 
lowest momentum bins. The thermal leakage injec- 
tion model emulates the process by which suprathermal 
particles well into the tail of the postshock Maxwellian 
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function of time. The Sedov Taylor similarity solutions for the outer shock, Rst and J7st are also shown for comparison. 
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distribution leak upstream across a quasi-parallel shock 
(Malkov & Volk 1998; Malkov & Drury 2001). This fil- 
tering process is implemented numerically by adopting 
a "transparency function", r osc (es,u), that expresses 
the probability of supra-thermal particles at a given 
velocity, v, leaking upstream through the postshock 
MHD waves (Kang et al. 2002). One free parameter 
controls this function; namely, e# = Bq/B±, which 
is the inverse ratio of the amplitude of the postshock 
MHD wave turbulence B± to the general magnetic field 
aligned with the shock normal, Bo (Malkov & Volk 
1998). Plasma hybrid simulations and theory both sug- 
gest that 0.25 < eb < 0.35, so that the model is well 
constrained. However, such large values of €b lead to 
very efficient initial injection and most of the shock en- 
ergy is quickly transferred to the CR component for 
strong shocks considered here (50 < M s < 300 ini- 
tially), causing a numerical problem at the very early 
stage of simulations. So we adopted smaller values, 



€b = 0.16 — 0.2 in this study. Dependence of the CR 
injection and acceleration on this parameter will be dis- 
cussed below. 

(d) A Bohm-like Diffusion Model 

Self-excitation of Alfven waves by the CR streaming 
instability in the upstream region is an integral part of 
the DAS at SNRs (Bell 1987; Volk et al. 1988; Lucek 
& Bell 2001). The particles are resonantly scattered by 
those waves, diffuse across the shock, and get injected 
into the Fermi first-order process. These complex in- 
teractions are represented by the diffusion coefficient, 
which is expressed in terms of a mean scattering length, 
A, and the particle speed, v, as n{x,p) — Xv/3. The 
Bohm diffusion model is commonly used to represent 
a saturated wave spectrum (i.e., A = r g , where r g is 
the gyro-radius), kb(p) = n n p 2 / (p 2 + l) 1 / 2 . Here 
Kn = mc 3 /(3eB) = 3.13 x lO^cmV 1 ^- 1 , and B^ 
is the magnetic field strength in units of microgauss. 
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Table 1. Model Parameters 



Model 


«ISM 


E 




es 


r a 


to 


u 


Po 




(cm -3 ) 


(10 51 ergs) 


iiG 




(PC) 


(years) 


(10 4 km s" 1 ) 


(10 6 crg cm 3 ) 


SI 


0.3 


1 


30 


0.16 


3.19 


255. 


1.22 


1.05 


S2 


0.3 


4 


30 


0.16 


3.19 


127. 


2.45 


4.20 


S3 


0.3 


1 


5 


0.16 


3.19 


255 


1.22 


1.05 


S4 


0.003 


1 


30 


0.16 


14.8 


1182. 


1.22 


1.05 (-2) 


S5 


0.3 


1 


5 


0.2 


3.19 


255 


1.22 


1.05 



Because of the steep momentum dependence for non- 
relativistic particles (p <C 1), simulations with a Bohm 
diffusion model require extremely fine grid resolution 
around the shock where freshly injected CRs are con- 
centrated. Instead we adopt a Bohm-like diffusion co- 
efficient that includes a weaker non-relativistic momen- 
tum dependence, 

K(r,p) = K n -p(j^J . (7) 

Previous studies showed that simulations using these 
two types of diffusion coefficient produced very simi- 
lar results (Kang et al. 2001). The assumed density 
dependence for n accounts for compression of the per- 
pendicular component of the wave magnetic field and 
also inhibits the acoustic instability in the precursor 
of highly modified CR shocks (Kang, Jones, & Ryu, 
1992). Hereafter we use the subscripts '0', '1', and '2' 
to denote conditions far upstream of the shock, imme- 
diately upstream of the gas subshock and immediately 
downstream of the subshock, respectively. 

III. Simulations of Sedov- Taylor Blast Waves 

For a supernova remnant propagating into a uniform 
ISM, the CR acceleration takes place mostly during 
free expansion and Sedov- Taylor (ST hereafter) stages, 
since the shock slows down significantly afterward. Fig. 
1 shows the evolution of a typical SNR calculated by 
a hydrodynamics code without the CR pressure terms. 
This demonstrates that the ST solution is established 
only after the inner reverse shock is reflected at the 
center at t/t Q <~ 7 (see below for the definition of nor- 
malization constants). Before that time, the reverse 
shock is strong and dynamically important. In our 
simulations, however, we will ignore the reverse shock, 
because the current version of CRASH code can treat 
only one shock. Application of our AMR algorithm for 
multiple spherical shocks is not simple, since it requires 
multiple, comoving grids. The CR acceleration at the 
reverse shock is thought to be not important, because 
the kinetic energy passed through the reverse shock is 
relatively small. Also adiabatic losses by CRs acceler- 
ated early on in the interior and then advected outward 
through the ST phase would generally be very large. 

On the other hand, Fig. 1 shows that the evolution 



of the outer shock speed, U ou t(t) can be approximated 
by the Sedov solution f7 ST oc (t/t )- 2 / 5 for t/t a > 0.2. 
In order to take account for the CR acceleration from 
free expansion stage through ST stage, we begin the 
calculations with the ST similarity solution at t/t Q = 
0.5. In principle wc could start the simulations from 
t/t a <~ 0.1. Such simulations, however, require very 
long computational time, since we need to include the 
extremely hot region with fast sound speeds near the 
explosion center. We carried out one model from t/t Q = 
0.2 and compared it with the case started from t/t Q = 
0.5 in the next section. Since the total CR energy gain 
is proportional to the kinetic energy passed through 
the shock, E sw — 2ir J p U^ ut Rl ut dt, we expect our 
calculations should capture the key aspects of the CR 
acceleration at the outer SNR shock. 

We terminated the simulation at t/t a — 10, while 
the ST stage ends typically when the shock slows down 
to Ust < 300 km s -1 around t/t Q <~ 3000. However, it 
has been shown that the highest momentum, p max , is 
achieved by the end of the free expansion stage and the 
transfer of explosion energy to the CR component oc- 
curs mostly during the early ST stage (e.g., Bcrczhko 
et al. 1997). We will show in the next section that 
the shock properties and the CR acceleration efficiency 
reach roughly time asymptotic values for t/t a > 1. Also 
it was suggested that non-linear wave damping and the 
wave dissipation due to ion-neutral collisions may sup- 
press the MHD waves significantly at the late ST stage, 
leading to fast particle diffusion and inefficient acceler- 
ation. (V61k et al. 1998; Ptuskin & Zirakashvili 2003) 

IV. Results 

(a) SNR Model Parameters 

We consider a Type la supernova explosion with the 
ejecta mass, M e j = 1.4Mq, in a warm or hot ISM with 
a uniform density. Model parameters are summarized 
in Table 1. The fiducial model, labeled SI in Table 1, 
has the explosion energy, E a = 10 51 ergs, and the back- 
ground density, riisM = 0.3 cm~ 3 . The pressure of the 
background gas is set to be Pism ~ 10~ 12 erg cm~ 3 , 
which determines the sound speed of the upstream gas 
and so the Mach number of the SNR shock. Recent 
X-ray observations of young SNRs indicate a magnetic 
field strength much greater than the mean ISM field of 
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Fig. 2. — Evolution of model SI SNR expanding into a uniform ISM at t/t a = 0.5, 1., 6., and 10. The model parameters 
are M ej = 1.4Mq, E a = 10 51 ergs, wism = 0.3cm~ 3 , and £? M = 30. The injection parameter for thermal leakage injection, 
es = 0.2. The lower left panel shows the volume integrated CR spectrum, G(p) = J f(r,p)p 4 r 2 dr. The initial condition at 
t/t = 0.5 (solid line) is set by the Sedov- Taylor similarity solution. 



5/iG, or values expected by compression of that field 
{e.g., Berezhko et al. 2003; Volk et al. 2005). It is be- 
lieved that the magnetic field upstream from the shock 
is amplified by the CR streaming instability in the pre- 
cursor region (Bell 1978; Lucek & Bell 2001). Thus, to 
represent this effect we take B — 30/iG as the fiducial 
field strength. The strength of magnetic field deter- 
mines the size of diffusion coefficient, n n , and the drift 
speed of Alfven waves relative to the bulk flow. The 
Alfven speed is given by va — va,o(p/ Pa)^ 1 ^ 2 where 
va,o = (1.8 kms -1 )i? At /^/nisM- The second model, 
S2, assumes a higher explosion energy, while the third 
model, S3, assumes the ISM magnetic field, rather 
than the amplified field. Model S4 assumes a hot ISM 
(T w 10 6 K), while all other models assume a warm 
ISM (T fa 10 4 K). For models S1-S4 e B = 0.16 is 
adopted. Model S5 has the same parameters as model 
S3 except es — 0.2, which allows a higher injection 



rate. 

The physical quantities are normalized, both in the 
numerical code and in the plots below, by the following 
constants: 



(3M ej 


\ 1/3 








\ 1/2 


\ E 




u = 


r /t c 



p = (2.34 x 10 - 24 gcm- 3 )ni SM , 

Po = PoU 2 . 

These values are also given in Table 1 for reference. 
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Fig. 3. — Immediate pre-subshock density, p\, post-subshock density, p2, shock Mach number, M B , post-subshock CR and 
gas pressure in units of the ram pressure of the unmodified Sedov- Taylor solution, poUg T oc (t/t )~ 6 ^ 5 , and the CR injection 
parameter, £, are plotted for models SI and S2. 



(b) Remnant Evolution 

Fig. 2 shows the evolution of model SI for t/t = 
0.5 — 10. The solid lines are for the initial condition 
which is assumed be the Sedov Taylor similarity solu- 
tion (R ST /r = Ut/to) - 6 , U ST /u = 0.6Ut/t o )-°-\ 
where £ s = 1.15167) extrapolated to t/t Q — 0.5. Ini- 
tially there is no pre-existing CRs and so all CR par- 
ticles are freshly injected at the shock. The CR pres- 
sure becomes dominated over the gas pressure and the 
density compression across the total shock becomes 
P2/P0 ~ 7 after t/t Q ~ 1. As the shock slows down and 
the CR pressure increases, the subshock Mach number 
decreases and the thermal population cools down, re- 
sulting in less efficient particle injection at low energies. 
We note that the inner boundary of the simulation grid 
moves out with the expanding comoving grid. 

We repeated the same calculation starting from an 
earlier time, t/t Q — 0.2, to explore how the CR ac- 



celeration during the early free expansion stage would 
affect the results at late ST stage. The total CR en- 
ergy accelerated by t/t Q = 10 differs about 3 % and 
the CR spectrum extends to slightly higher p max in the 
simulation started from t/t = 0.2, as one would antic- 
ipate from the longer acceleration interval. The small 
differences indicate that omitting the evolution before 
t/t a = 0.5 does not affect the overall conclusions of this 
work. 

(c) CR Injection and Acceleration 

The efficiency of the particle injection is quantified 
by the fraction of particles swept through the shock 
that have been injected into the CR distribution: 

= 1 47rr 2 dr / 47r/ CR (p, r, f)p 2 dp 
J 47rr 2 nou s dt 

where /or is the CR distribution function, no is the 
particle number density far upstream and r s is the 
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Fig. 4. — Pre-subshock density, pi, post-subshock density, p2, shock Mach number, M s , post-subshock CR and gas pressure 
in units of the ram pressure of the unmodified Sedov- Taylor solution, polIg T oc (t/t )~ 6 ^ 5 , and the CR injection parameter, 
£, are plotted for models S3 and S4. The left panels also show model S5 with es = 0.2 (long dashed lines and dotted lines). 



shock radius. 

Figs. 3-4 show the evolution of shock properties such 
as the compression factors, subshock Mach number, 
postshock pressures, and the injection fraction for all 
five models. Most of these quantities approach to time- 
asymptotic values. The shock Mach number is the key 
parameter that determines the CR injection and ac- 
celeration efficiency. Model S4 has much lower initial 
Mach number, so it shows less injection rate (£ ~ 10 -5 ) 
compared to the other models (£ J> 10~ 4 ). The com- 
pression factor in the precursor, pi/pa ~ 2 — 3, while 
the compression factor across the total shock structure 
varies somewhat, P2/P0 ~ 7.0, 8.2, 10., and 4.8 for 
models S1-S4, respectively. The postshock CR pres- 
sure relative to the ram pressure of the Sedov solution 
is P c , 2 /pof/| T ~ 0.35, 0.39, 0.42, and 0.25 for models 
S1-S4, respectively. 

We note that model S4 with B = 5/iG has the 
higher CR pressure than the fiducial model SI with 



B = 30/iG. While the mode with higher B has smaller 
diffusion coefficient and so the CR spectrum extends to 
higher p max (see Fig. 5), the reduction of the injection 
and acceleration due to the wave drift leads to less CR 
energy. 

Model S5 has a larger value of cb = 0.2 than model 
S3 does, but otherwise they have the same model pa- 
rameters. Weaker turbulence (larger e_e) lead to higher 
thermal leakage, so the injection rate £ is about 15 
times larger than that of models S1-S3. This in turn 
results in slightly higher CR pressure and more signifi- 
cant precursor compression. However, it demonstrates 
that the CR acceleration efficiency depends only weakly 
on the injection rate. 

Fig. 5 shows the volume integrated CR spectrum, 
G(p) = f 4irg(p)r 2 dr in code units and the slope, 
q = —d(bxGp)/d\a.p — 2. The maximum momentum 
depends on the SNR parameters as p max oc u s t/K n oc 
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Fig. 5. — Integrated CR number, G p = J g(p)r 2 dr, in arbitrary units and the slope, q = — d(ln G p )/d hip — 2, are shown at 
£/to = 1, 3, 6, and 10. See Table 1 for model parameters for models Sl-4. 



b h e oPism ■ Tne values of log(p max ) w 6.5, 6.8, 5.7, 
and 7.0 at </t D = 10 for models S1-S4, respectively, are 
roughly consistent with this relation. 

The CR spectra for models S1-S3 with high initial 
Mach numbers show the canonical nonlinear concave 
curvature for p > 1, which comes from the large den- 
sity compression factor (pz/po ^> 4) due to the domi- 
nant CR pressure. For these models the slope q is close 
to 2 at low energies, but flattens to ~ 1.6 at high ener- 
gies below the upper momentum cutoff. For model S4 
with lower initial Mach number the density compres- 
sion is only P2/P0 ~ 4.8, so the flattening of G(p) at 
high momenta is not so prominent. Instead, the de- 
crease of the CR injection rate due to the weakening 
subshock reduces the particle numbers at low momen- 
tum, leading to actually flatter spectra there. 

Fig. 6 shows the integrated energies, Ei/E a = 
4-7T J eir 2 dr, where eth, ekin, and ec are the density 
of thermal, kinetic and cosmic ray energy, respectively. 
The kinetic energy reduces only slightly and is similar 



for all models. The total CR energy accelerated up to 
tjt = 10 is E c /E = 0.55,0.57,0.62 and 0.39 for mod- 
els S1-S4, respectively. So models S2 (larger E a ) and 
S3 (smaller B^) are a bit more efficient in transferring 
the SN explosion energy to the CR energy, compared 
to the fiducial model. As mentioned earlier, compari- 
son between models SI and S3 indicates that strongly 
amplified magnetic field carries faster Alfven waves and 
leads to less efficient CR acceleration, even though the 
particles are accelerated to higher energies. Also com- 
parison between models S3 and S5 shows that the total 
CR energy increases only slightly, less than 3 %, even 
though the injection rate increases by a factor of about 
15 with a larger value of es for model S5. In model S4 
with a hotter ISM, the shock is weaker and so the CR 
injection and acceleration are less efficient. 

Our simulations imply that on average about 10 -4 — 
10~ 3 of the incoming particles are injected to the CR 
population at the shock front and up to 50 % of the 
SN explosion energy can be transferred to CRs during 
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Fig. 6. — Integrated thermal, kinetic and CR energies inside the simulation volume as a function of time. See Table 1 for 
model parameters for models S1-S5. The long dashed line in the lower left panel shows the volume integrated CR energy 
for model S5. 



the ST stage, if the magnetic field direction is radial 
(i.e., quasi-parallel field). This last result is consis- 
tent with the calculations previously done by Berezhko 
and collaborators using a different numerical scheme 
(e.g., Berezhko & Volk 1997). The CR injection rate, 
however, probably depends strongly on the angle be- 
tween the magnetic field and the shock normal direc- 
tion. In a more realistic magnetic field geometry, where 
a uniform ISM field is swept by the spherical shock, 
only 10-20 % of the shock surface has a quasi-parallel 
field geometry (Volk et al. 2003). In the shock surface 
region where the field is perpendicular, the injection 
rate is expected to be reduced and so the CR accel- 
eration efficiency would be smaller. Thus the CR en- 
ergy conversion factor averaged over the entire shock 
surface could be significantly smaller than 50 %, per- 
haps about 10 %. On the other hand, Giacolne (2005) 
showed that the protons can be injected efficiently even 
at perpendicular shocks in fully turbulent fields due to 
field line meandering. In such case the injection rate at 



perpendicular shocks may not be significantly 
compared to parallel shocks. 



smaller, 



V. SUMMARY 



The evolution of cosmic ray modified shocks de- 
pends on complex interactions between the particles, 
waves in the magnetic field, and underlying plasma 
flow. We have developed numerical tools that can em- 
ulate some of those interactions and incorporated them 
into a standard numerical scheme for gasdynamic prob- 
lems. Specifically, diffusive shock acceleration can be 
followed by a kinetic approach (i.e., CRASH code) in 
which a diffusion convection equation for the CR distri- 
bution function is solved with an appropriate diffusion 
coefficient (Kang & Jones 1991; Kang et al. 2001). The 
injection of CR particles can be treated by a thermal 
leakage injection model with a transparency function 
T e sc(cB,v) (Kang et al. 2002). Drifts of resonantly- 
scattering Alfven waves and heating of the thermal 
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plasma due to the dissipation of those waves in the 
precursor region has been included through a simple 
model (Jones 1993; Kang & Jones 2006). 

In the present paper, we applied the spherical CRASH 
code to the problem of the CR acceleration at the rem- 
nant shock from typical Type la SNe propagating into 
a uniform interstellar medium. The main results of our 
simulations can be summarized as follows: 

1) Since the CR injection and acceleration depend 
primarily on the shock Mach number, the temperature 
of the ISM is an important factor. SNRs in a warm 
ISM inject about 10~ 4 — 10 -3 of the particles passed 
through the shock and transfer about 50% of the ex- 
plosion energy to the CR component. SNRs in a hot 
ISM inject 10 times less particles, but they still transfer 
about 40 % of the explosion energy to CRs. 

2) At sources, the hadronic CR spectrum should 
be consistent with a power-law function, N(E)dE oc 
E~ 2 dE up to p/m p c <~ 10 2 , but it may gradually flat- 
ten to E~ 16 toward the Knee energy and beyond up 
to E ~ 10 16 Z eV (p/m p c ~ 10 7 ). 

3) About 40-50 % of the explosion energy can be 
transfered to CRs, if the magnetic field is radial (i.e., quasi- 
parallel shocks). However, more realistic consideration 

of the field geometry relative to the spherical shock sur- 
face could lead to much smaller energy conversion rate. 
So a conservative estimate could be order of 10 %. 

4) If the magnetic field is amplified in the precursor 
region due to the streaming instability, as indicated by 
recent X-ray observations of young SNRs, the particles 
are accelerated to higher energies due to smaller dif- 
fusion coefficient. However, faster Alfven waves in the 
precursor tend to reduce the CR injection and acceler- 
ation efficiencies. 

In conclusion, the galactic cosmic rays possibly up 
to 10 16 Z eV could be originated from supernova rem- 
nants. This result is quite robust, regardless of SNR 
model parameters and details of microphysics involved 
in the injection process, provided that the Bohm diffu- 
sion is valid at young SNRs of several thousand years 
old. 

ACKNOWLEDGEMENTS 

This work was supported for two years by Pusan 
National University Research Grant. The author would 
like to thank T.W. Jones for helpful comments on the 
paper. 

REFERENCES 

Aharonian et al. , H.E.S.S. collaboration, 2004, High-energy 
particle acceleration in the shell of a supernova remnant, 
Nature, 432, 75 

Bamba, A., Yamazaki, R, Ueno, M. & Koyama, K., 2003, 
Small-Scale Structure of the SN 1006 Shock with Chan- 
dra Observations, ApJ, 589, 827 

Bell, A. R., 1978, The acceleration of cosmic rays in shock 
fronts. I, MNRAS, 182, 147 



Berezhko, E. G., Ksenofontov, L. T., & Yelshin, V., 1995, 
Efficiency of CR acceleration in supernova remnants, Nu- 
clear Physic B., 39, 171. 

Berezhko, E. G, & Volk, H. J., 1997, Kinetic theory of 
cosmic rays and gamma rays in supernova remnants. I. 
Uniform interstellar medium, Astropart. Phys. 7, 183 

Berezhko, E. G, Ksenofontov, L. T., & Volk, H. J., 2003, 
Confirmation of strong magnetic field amplification and 
nuclear cosmic ray acceleration in SN 1006, A&Ap, 423, 
Lll 

Berezhko, E. G, & Volk, H. J., 2006, Theory of cosmic ray 
production in the supernova remnant RX J1713. 7-3946, 
A&Ap, 451, 981 

Blandford, R. D., & Ostriker, J. P., 1978, Particle acceler- 
ation by astrophysical shocks, ApJ, 221, L29 

Blandford, R. D., & Eichler, D., 1987, Particle acceleration 
at astrophysical shocks - a theory of cosmic-ray origin, 
Phys. Rept., 154, 1 

Drury, L. O'C, 1983, An Introduction to the Theory of 
Shock Acceleration of Energetic Particles in Tenuous 
Plasmas, Rept. Prog. Phys., 46, 973 

Drury, L. O'C, Ellison, D. E., Aharonian, F. A. et al. , 
2001, Test of galactic cosmic-ray source models - Working 
Group Report, Space Science Reviews, 99, 329 

Giacalone, J., 2005, The Efficient Acceleration of Thermal 
Protons by Perpendicular Shocks ApJ, 628, L37 

Jones, T. W., 1993, Alfven wave transport effects in the 
time evolution of parallel cosmic-ray-modified shocks 
ApJ, 413, 619 

Kang, H., & Jones, T. W., 1991, Numerical studies of dif- 
fusive particle acceleration in supernova remnants, MN- 
RAS, 249, 439 

Kang, H., Jones, T. W., & Ryu, D., 1992, Acoustic insta- 
bility in cosmic ray mediated shocks, ApJ, 385, 193 

Kang, H., Jones, T. W., LeVeque, R. J., & Shyue, K. M., 
2001, Time Evolution of Cosmic-Ray Modified Plane 
Shocks, ApJ, 550, 737 

Kang, H., Jones, T. W., & Gieseler, U. D. J., 2002, Numer- 
ical Studies of Cosmic-Ray Injection and Acceleration, 
ApJ, 579, 337 

Kang, H., & Jones, T. W., 2002, Acceleration of Cosmic 
Rays at Large Scale Cosmic Shocks in the Universe, Jour- 
nal of Korean Astronomical Society, 35, 159 

Kang, H., 2003, Cosmic Ray Acceleration at Cosmological 
Shocks: Numerical Simulations of CR Modified Plane- 
Parallel Shocks Journal of Korean Astronomical Society, 
36, 1 

Kang, H., & Jones, T. W., 2006, Numerical studies of dif- 
fusive shock acceleration at spherical shocks Astropart. 
Phys. 25, 246 

Koyama, K, Petre, R., Gotthelf, E. V. et al. , 1995, Ev- 
idence for Shock Acceleration of High-Energy Electrons 
in the Supernova Remnant SN:1006, Nature, 378, 255 

Lagage, P. O., & Cesarsky, C. J., 1983, The maximum 
energy of cosmic rays accelerated by supernova shocks, 
A&Ap, 118, 223 



Cosmic Ray Acceleration at SNRs 



Lucek, S. G., & Bell, A. R., 2000, Non-linear amplifica- 
tion of a magnetic field driven by cosmic ray streaming, 
MNRAS, 314, 65 

Malkov, M. A., & Drury, L. O'C, 2001, Nonlinear theory 
of diffusive acceleration of particles by shock waves, Rep. 
Progr. Phys. 64, 429 

Malkov, M. A., & V61k, H. J., 1998, Diffusive ion acceler- 
ation at shocks: the problem of injection, Adv. Space 
Res. 21, 551 

Ptuskin, V. S., & Zirakashvili, V. N., 2003, Limits on dif- 
fusive shock acceleration in supernova remnants in the 
presence of cosmic-ray streaming instability and wave 
dissipation, A&Ap, 403, 1 

Skilling, J., 1975, Cosmic ray streaming. I - Effect of Alfven 
waves on particles, MNRAS, 172, 557 

Volk, H. J., Berezhko, E. G., & Ksenofontov, L. T., 
2003, Variation of cosmic ray injection across supernova 
shocks, A&Ap, 409, 563 

Volk, H. J., Berezhko, E. G., & Ksenofontov, L. T., 2005, 
Magnetic field amplification in Tycho and other shell- 
type supernova remnants A&Ap, 433, 229 

Volk, H. J., Zank, L. A., & Zank, G. P., 1988, Cosmic 
ray spectrum produced by supernova remnants with an 
upper limit on wave dissipation, A&Ap, 198, 274 



